EMPIRICAL  ANALYSIS  OL  A  3-D  PAKISTAN  REGIONALIZED  VELOCITY  MODEL 
FOR  IMPROVED  SEISMIC  EVENT  LOCATIONS 


Michelle  Bernard-Johnson,  Shirley  Rieven,  Carolynn  Vincent,  Delaine  Reiter, 

Weston  Geophysical 

William  Rodi, 

Earth  Resources  Laboratory,  Massachusetts  Institute  of  Technology 

Sponsored  by  U.S.  Department  of  Defense 
Defense  Threat  Reduction  Agency 
Contract  No.  DSWA01-98-C-0143 


ABSTRACT 


The  objective  of  this  project  is  to  improve  regional  event  location  in  Pakistan  and  the  surrounding  regions  by 
developing  a  3-D  velocity  model  for  the  crust  and  upper  mantle.  The  CTBT  goal  of  1000  km2  location  accuracy 
for  small  events  is  dependent  upon  the  availability  of  such  models,  which  can  be  used  to  compute  accurate  travel 
times  of  regional  seismic  phases.  Our  approach  is  to  refine  an  a  priori  3-D  velocity  model  through  joint 
nonlinear  tomography  and  multiple  event  location  applied  to  regional  earthquake  arrival  time  data.  Previously, 
we  presented  our  a  priori  model,  which  is  constructed  on  a  l°-by-l°  grid  to  a  depth  of  approximately  100  km. 
The  preliminary  model  exhibits  complex  structural  variation  in  the  Pakistan  region  with  crustal  thickness  ranging 
from  25  km  to  75  km,  a  clear  deviation  from  1-D  global  models  such  as  the  IASP91  model  (Kennett  and  Engdahl, 
1991). 

We  have  generated  3-D  P,  Pg  and  Pn  travel  time  tables  at  regional  distances  by  applying  the  finite-difference 
raytracing  method  of  Podvin  and  Lecomte  ( 1991)  to  our  a  priori  model.  To  test  our  model  empirically,  we  have 
relocated  a  suite  of  44  Calibration  Event  Bulletin  events  using  P  and  Pn  arrival  times  from  14  stations  in  the 
Pakistan  region,  reported  to  the  International  Seismological  Centre  (ISC).  Locations  obtained  with  the  IASP91 
travel  time  tables  were  compared  to  locations  obtained  with  our  3-D  travel  time  tables.  While  our  3-D  model 
results  in  smaller  mis  residuals  for  29  of  44  events,  the  tests  are  not  definitive  since  ground-truth  locations  are  not 
available  and  the  depths  of  the  events  are  poorly  constrained  by  the  regional  data. 

We  are  currently  expanding  our  database  to  include  55,000  phase  arrivals  from  over  4000  events  recorded  at  77 
stations  in  the  study  region,  which  we  are  expanding  to  include  additional  parts  of  southwest  India.  The  arrival 
time  picks  at  these  and  teleseismic  stations  were  reassociated,  and  the  events  relocated,  by  Engdahl  el  al.  (1998), 
leading  to  high  quality  location  estimates.  Approximately  ten  percent  of  the  event  locations  are  classified  as  GT5 
(5  km  location  accuracy  )  or  better.  Using  the  reassociated  P  times  at  the  regional  stations,  we  intend  to  relocate  a 
subset  of  the  events  using  both  the  IASP91  and  our  3-D  travel  time  tables.  A  more  extensive  evaluation  of  our  3— 
D  model  will  then  be  possible  by  comparing  rms  residuals,  as  well  as  comparing  our  locations  to  those  obtained 
by  Engdahl  et  al.  Furthermore,  the  residuals  resulting  with  our  3-D  tables  will  comprise  the  primary  data  set  for 
tomographic  updating  of  our  a  priori  model.  We  are  analyzing  the  spatial  pattern  of  these  residuals  and  raypath 
coverage  in  anticipation  of  performing  the  3-D  tomography. 
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OBJECTIVE 


Precise  location  of  seismic  events  in  the  context  of  the  Comprehensive  Nuclear-Test-Ban  Treaty  (CTBT) 
requires  accurate  velocity  models.  The  purpose  of  this  research  is  to  construct  a  regional  velocity  model  of  the 
crust  and  upper  mantle  structure  to  improve  seismic  event  locations  for  the  area  between  25  and  40  degrees  North 
and  55  and  80  degrees  East.  This  area  lies  between  the  eastern  border  of  Iran  and  the  western  edge  of  the  Tibetan 
Plateau.  To  date  we  have  developed  a  database  of  regional  seismic  events  and  supplemental  information  that  will 
be  used  to  construct  and  refine  the  velocity  model  through  the  application  of  tomographic  techniques.  We 
generate  travel  times  through  the  region  using  a  finite  difference  travel  time  modeling  algorithm  (Podvin  and 
Lecomte  ,  1991 )  and  use  our  grid  search  algorithm  for  3-D  hypocenter  relocation.  We  are  developing  a  joint 
inversion  technique  that  uses  multiple  event  relocation  (incorporating  input  from  our  database  of  phase  arrivals)  to 
constrain  iterative  perturbations  of  the  velocity  model.  This  joint  inversion  will  provide  the  “best”  fit,  in  the 
least-squares  sense,  between  the  observed  and  predicted  arrival  times.  The  end  result  will  be  a  refined  3-D 
velocity  model  and  improved  multiple  hypocenter  estimates.  In  this  report  we  describe  the  progress  we  have 
made  towards  building  a  high  quality  regional  model  during  the  second  year  of  this  three  year  project. 

RESEARCH  ACCOMPLISHED 


Phase  Arrival  Database  and  3-D  Velocity  Model  Development  and  Testing 

Two  primary  objectives  during  the  first  two  years  of  this  project  have  been  to  collect  and  refine  a  high  quality 
database  of  phase  arrivals  and  to  develop  a  good  3-D  preliminary  velocity  model  over  the  region.  Both  the 
database  and  the  3-D  velocity  model  can  be  used  in  tomographic  studies  as  well  as  other  related  projects.  During 
this  reporting  period  we  have  addressed  gaps  in  the  ray-coverage  over  Pakistan  by  expanding  our  study  area 
boundaries  to  the  West  (covering  areas  of  Iran),  North  (including  areas  of  Uzbekistan,  Kyrgyzstan  and 
Kazahkstan)  and  East  (China  and  portions  of  northern  India)  where  additional  data  were  readily  available.  We  are 
also  expanding  our  region  to  the  South  to  include  areas  of  northwestern  India.  Below  we  detail  some  of  the 
development  and  testing  of  the  phase  arrival  database  and  3-D  velocity  model  that  has  occurred  during  the 
previous  year. 

Travel  Time  Database  Improvements 

Our  initial  effort  to  establish  a  database  of  regional  phase  arrivals  began  with  the  acquisition  of  bulletin  data  for 
approximately  700  seismic  events.  These  phase  arrivals  were  obtained  from  three  online  sources,  including  the 
Prototype  International  Data  Center  (pIDC),  the  International  Seismological  Centre  (ISC),  and  the  Incorporated 
Research  Institutes  for  Seismology  (IRIS).  During  the  early  stages  of  our  research,  we  used  these  data  to  identify 
a  set  of  test  events  that  would  contribute  the  "observables"  portion  of  our  input  data  and  also  provide  a  set  of 
validation  events.  Because  these  events  have  been  located  by  multiple  organizations,  the  redundancy  of 
hypocenter  estimates  and  the  availability  of  depth  phase  data  provided  a  measure  of  the  accuracy  of  the  combined 
event  catalog.  Upon  comparison  of  different  estimates,  it  became  apparent  that  the  relative  quality  of  the  location 
estimates  exhibited  significant  variation,  particularly  in  depth.  The  pIDC  reported  very  few  ground  truth  (GT) 
events  with  location  accuracies  better  than  25km  (GT25).  This  lack  of  validation  events  presents  a  significant 
challenge  in  establishing  a  demonstrable  advantage  in  using  a  3-D  model  over  the  standard  1-D  models  currently 
used.  In  this  performance  period,  we  have  investigated  methods  of  addressing  these  concerns,  and  give  details  of 
some  of  these  procedures  later  in  the  paper. 

We  have  expanded  our  in-house  database  of  phase  arrivals  by  incorporating  the  results  of  very  recent  hypocenter 
relocation  studies.  In  particular,  we  compared  the  relative  hypocenter  accuracies  of  the  ISC  catalog  with  those 
reported  by  Engdahl  el  al.  (1998).  In  that  study,  the  authors  reassociated  the  ISC  arrivals  and  then  relocated  the 
events  using  the  AK135  velocity  model.  Our  original  database  reflected  the  hypocentral  estimates  as  reported 
only  by  the  three  locating  agencies  previously  noted.  For  instance,  the  pIDC’s  Ground  Truth  Database  for  the 
Pakistan  region  identified  only  60  GT  events  during  its  period  of  data  availability.  Of  those,  only  7  are  reported  as 
being  located  with  hypocenter  accuracies  better  than  GT  25  (one  event,  the  1974  Indian  nuclear  test,  is  listed  as  a 
GT0,  while  six  events  are  classified  as  GT1,  including  the  1998  India  and  Pakistan  nuclear  tests).  By  comparison, 
the  Engdahl  el  al.  (1998)  results  report  2,181  events  qualifying  as  GT10  or  better  (including  8  explosions),  of 
which  there  are  426  events  qualifying  as  GT  5  or  better.  The  availability  of  these  additional  well  located  events 
improves  our  ability  to  verify  the  quality  of  our  tomographically  derived  velocity  model.  In  addition,  the  Engdahl 


et  al.  (1998)  database  offers  an  additional  advantage,  because  P  waves  were  reclassified  as  Pn  when  appropriate. 
This  provides  an  important  data  distinction  that  we  take  advantage  of,  since  in  our  method  we  invert  for  these 
secondary  arrivals  using  independently  derived  travel  time  calculations. 

In  summary,  our  database  has  expanded  from  7,599  phase  arrivals  to  54,770  P.  Pn,  and  Pg  arrivals.  These  events 
are  distributed  in  depth,  with  the  predominant  proportion  (76.5  percent)  being  located  within  the  top  100  km  of  the 
Earth's  surface.  The  magnitudes  of  these  events  range  from  a  high  of  6.5  (mb)  to  a  minimum  of  3.0.  More  than 
80  percent  of  these  events  have  at  least  10  regional  stations  reporting  phase  arrivals  A  total  of  77  regional 
stations  contributed  arrivals  to  the  expanded  database.  We  are  currently  in  the  process  of  evaluating  the  azimuthal 
gaps  for  these  events.  We  will  continue  our  effort  to  identify,  evaluate  and  incorporate  phase  data  from  events 
with  high  quality  location  estimates  .  To  ensure  an  independent  test  of  our  tomographic  model,  we  will  establish 
a  set  of  our  highest  quality  events  to  be  used  solely  for  verification  purposes..  The  remainder  will  be  used  in  the 
tomographic  inversion. 

Velocity  Model  Improvements 

High  density,  well  distributed  ray  coverage  is  crucial  to  the  success  of  a  tomographic  inversion,  because  it 
constrains  the  possible  velocities  for  large  portions  of  the  model  grid.  We  have  increased  the  size  of  our  study 
area  to  incorporate  rays  that  travel  significant  distances  within  our  contractual  region  but  are  not  contained 
entirely  within  it.  This  increase  in  sudy  size  has  required  the  collection  of  additional  geologic  and  geophysical 
information  to  expand  the  model  to  the  South  into  sections  of  India.  There  are  a  number  of  stations  in  this  region 
that  will  provide  travel  times  through  the  southeastern  quadrant  where  the  southwestern  portion  of  the  Indus 
Suture  defines  the  convergence  of  the  Indian  plate  with  the  Eurasian  plate.  Incorporating  rays  that  cross  this 
region  is  especially  important  for  the  purposes  of  CTBT  monitoring  because  seismic  events  that  occur  near  test 
sites  would  propagate  rays  through  this  region.  Accurate  hypocenter  location  estimates  would  rely  heavily  on  the 
quality  of  the  velocity  model  available  for  this  region.  In  the  following  sections  we  cover  some  of  the  procedures 
we  have  implemented  to  test  the  validity  of  our  input  3-D  regional  velocity  model. 

Velocity  Model  Testing  Using  Reciprocity 

On  14  Feb  1977,  a  magnitude  5.2  earthquake  (Table  1)  occurred  in  the  region  near  Nilore,  Pakistan.  The  event 
forms  an  interesting  test  of  reciprocity  for  the  velocity  model  near  station  NIL.  We  collected  arrival  time  data 
from  24  stations  within  our  model  region  that  recorded  the  event,  and  then  compared  the  observed  travel  time 
residuals  at  these  stations,  with  respect  to  the  Jeffreys-Bullen  ( JB )  global  model,  to  the  predicted  travel  time 
residuals  computed  using  our  preliminary  3-D  model  (Figure  1 ).  The  results  show  that  our  3-D  model  predicts 
the  regionally  observed  travel  times  better  than  the  JB  model  by  a  factor  of  1.5.  We  expect  that  the  largest 
residuals  between  the  predicted  and  observed  travel  times  (Figure  2)  will  decrease  following  a  full  tomographic 
inversion  of  our  dataset. 


Date 

OT 

Lot 

Long 

Depth  (km) 

Mag 

14  Feb  1977 

00:22:38 

33.5967 

73.2669 

27 

5.2 

Table  1.  ISC  Origin  data  for  an  earthquake  that  occurred  beneath  station  NIL.  The  ISC  listed  the  depth  as  27 
km;  however,  detailed  aftershock  studies  in  the  epicentral  region  suggest  the  14  Feb  1977  main  shock  occurred  at 
a  depth  of  12-20  km  (Seeber  and  Armbruster,  1979).  To  accommodate  this  depth  range,  we  computed  the  SSSC 
for  station  NIL  at  15  km. 
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Figure  1.  Observed  Pg  and  P n  travel  times  at  24  regional  stations  that  recorded  the  14  Feb  1977  event  underneath  NIL.  The 
predicted  travel  times  were  calculated  using  the  Weston  3-D  velocity  model  for  Pakistan  and  were  found  to  model  the 
observed  data  by  a  factor  of  1.5  times  better  than  the  fit  provided  by  the  Jeffreys-Bullen  model. 
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Figure  2  Difference  between  observed  and  predicted  travel  times  as  plotted  at  the  stations  that  recorded  the  Nilore  (NIL) 
earthquake. 

Velocity  Model  Testing  Using  Jackknifing 

We  are  interested  in  establishing  a  baseline  understanding  of  how  our  preliminary  input  model  performs  in 
comparison  to  the  1-D  IASP91  model.  Comparing  rms  travel  time  residuals  from  event  locations  using  the  two 
different  models  has  shown  inconsistent  results.  The  jackknife  method  of  estimating  error  (Efron  and  Gong,  1983) 
provides  another  comparison  of  location  performance  between  our  3-D  model  and  IASP9 1 .  In  this  method,  we 
relocate  each  event  multiple  times  (as  many  times  as  there  are  observed  phase  arrivals  in  the  inversion),  and  for 
each  relocation,  we  exclude  one  phase  arrival.  Inspecting  the  spatial  distribution  of  all  of  the  relocation  estimates 
provides  a  measure  of  how  each  data  point  affects  the  inversion.  Large  changes  in  the  hypocenter  introduced  by  a 
particular  datum  can  indicate  either  poor  quality  phase  picking  or  it  can  indicate  a  poor  model  along  that  ray  path. 
A  good  model  and  well  picked  data  would  produce  a  full  set  of  relocation  estimates  that  are  tightly  clustered. 

We  applied  the  jackknife  method  to  four  events  from  the  Engdahl  et  al.{  1998)  database.  These  4  events  are  a 


subset  of  200  that  we  used  for  a  pre-tomographic  analysis  test  case  discussed  later  in  this  paper.  The  results  of 
the  4-event  jackknife  study  are  shown  in  Figure  3.  Events  1  and  2  have  24  and  22  arrivals,  respectively,  and 
show  smaller  hypocenter  changes  relative  to  events  3  and  4  which  have  only  10  arrivals.  The  rms  residuals  for 
events  1  and  3  were  smaller  for  our  3-D  model  (0.923  s  and  1.196  s,  respectively)  than  for  the  IASP91  model 
(1.071  s  and  1.594  s  ).  For  event  2,  the  residuals  were  similar  (0.844  versus  0.883).  The  residual  for  event  4  is 
larger  in  the  3-D  model  case  (1.572  s)  than  in  the  IASP91  case  (1.414  s).  The  results  of  the  jackknife  tests 
indicate  that,  in  all  cases,  our  3-D  model  produces  smaller  relocation  scatter  than  does  the  IASP91  model  (See 
Table  2  for  scatter  values).  This  indicates  that  even  our  initial  3-D  model  produces  a  more  robust  fit  to  the 
observed  data  than  does  the  1-D  model,  and  a  3-D  model  improved  by  tomographic  inversion  should  perform 
even  better. 


EVENT 

ID 

3-D 

Epicenter 
Scatter  (KM) 

IASP91 

Epicenter 

Scatter 

(km) 

3-D 

Depth 
scatter 
(km  ) 

IASP91 

Depth 

scatter 

(km) 

Number 
of  Data 

1 

0.6750 

0.7375 

1.6083 

1.9542 

24 

2 

0.5636 

0.6227 

1.2772 

1.2864 

22 

3 

8.2600 

12.610 

27.810 

48.790 

10 

4 

4.3600 

4.7500 

5.7000 

7.9900 

10 

Table  2.  Results  of  the  jackknife  study  of  four  events  from  the  Engdahl  et  al.  (1998)  database.  The  values 
indicate  average  scatter. 


map  view 


latitude  vs.  depth 


longitude  vs.  depth 


^  35.9 
1  35.8 
C  S  35.7 

0}  m 

>  |  35.6 

•J  35.5 
35.4  _ 


0 

0 

50 

50 

E 

m 

3  100 

Q. 

□ 

£  100 
CL 

□ 

Q 

150 

I 

Q 

150 

200 

200 

70.4  70.6  70.8 

Longitude  (degrees) 


35.4  35.6  35.8 

Latitude  (degrees) 


70.4  70.6  70.8 

Longitude  (degrees) 


35.3 

|  35.2 

f  35.1 

■8  35 
3 

3  34.9 
34.8 


42.6 
1 42-5 

42.4 

jj  42.3 
1  42.2 
42.1 


0 

0 

50 

50 

E 

1 

□ 

3  100 

Q. 

□ 

£  100 
Q. 

□ 

Q 

a 

150 

« 

150 

% 

200 

200 

70.2  70.4  70.6 

Longitude  (degrees) 


34.8  35  35.2 

Latitude  (degrees) 

Or* — O - * - 


70.2  70.4  70.6 

Longitude  (degrees) 


72.8  73  73.2 

Longitude  (degrees) 


50 

E 

3  100 

Q. 

<D 

o 

150 

200  - 


0 

50 

r  ioo 

L 

) 

I 

150 


42.2  42.4  42.6 

Latitude  (degrees) 


73  73.2  73.4 

Longitude  (degrees) 


0 

□ 

0 

□ 

50 

50 

E 

3  100 

£  100 

a. 

C-T^X 

a 

a 

□  X 

150 

o  0 

150 

O  O 

>  m* 

xl*x 

200 

— . - ^ . - 

200 

- .  &  X  . - 

70.8  71  71.2 

Longitude  (degrees) 


Latitude  (degrees) 


70.8  71  71.2  71.4 

Longitude  (degrees) 


Figure  3.  Jackknife  test  results  for  four  events.  "X"s  represent  3-D  model  hypocenters,  "o"s  represent  1-D  IASP91  model 
hypocenters.  Singular  squares  represent  the  hypocenters  calculated  by  Engdahl  et  al.  (1998)  using  the  AK135  model. 


Preliminary  Tomographic  Analysis 


Figure  4  illustrates  the  five  main  components  in  our  tomographic  inversion  method  that  have  either  been 
completed  or  are  in  development.  The  "Condition  Model"  function  transforms  input  velocity  models  into  the 
proper  format  for  the  inversion.  Converting  from  spherical  to  Cartesian  coordinates  produces  a  station-centered 
sub-grid  which,  in  part,  reduces  errors  introduced  by  the  coordinate  transformation.  This  module 
reparameterizes  the  velocity  model  to  reduce  the  size  of  the  tomographic  inversion  and  constrain  the  inversion  to 
geologically  realistic  solutions.  The  "Predict  travel  time"  function  produces  travel  time  tables  for  each  station. 
Using  a  finite  difference  implementation  of  Huygen’s  principle  (Podvin  and  Lecomte,  1991;  Lomax,  1999),  we 
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Figure  4:  Schematic  overview  of  joint  tomographic  inversion  for  a  3-D  velocity  model  and  improved  seismic  event 
hypocenters.  The  degree  of  completion  of  a  particular  function  is  indicated  by  the  shading  of  the  block,  where  lighter  shades 
indicate  modules  which  require  further  development. 


Our  3-D  grid  search  algorithm  ("Invert  for  location" )  produces  hypocenters  for  each  event  in  the  tomographic 
inversion.  More  specifically,  it  obtains  the  maximum  likelihood  estimate  (MLE)  of  the  hypocentral  parameters 
and  the  probability  density  function  scale  factor  by  searching  a  grid  of  solution  hypocenters  through  subsets  of 
increasing  fineness.  The  final  MLE  of  the  hypocenter  is  taken  as  the  solution  achieving  the  largest  likelihood 
among  all  the  points  of  the  grids  searched.  We  used  this  algorithm  to  generate  the  estimates  shown  in  the 
jackknife  test  discussed  previously. 

In  order  to  use  the  Podvin  and  Lecomte  (1991)  travel  time  estimation  method  as  the  forward  modeling  core  of  the 
tomographic  inversion,  we  are  expanding  its  functionality  to  include  the  generation  of  accurate  spatial  coordinates 
of  the  rays  as  they  travel  from  the  source  to  the  receiver.  The  lack  of  accurate  ray  coordinates  is  a  deficiency  in 
the  original  version  of  the  original  that  must  be  remedied  prior  to  performing  tomographic  inversion.  The 
improved  ray  tracing  module  will  produce  ray  coordinates  for  each  input  station/event  phase  arrival  time,  which 
will  be  recalculated  after  every  update  of  the  velocity  model  until  the  fit  to  the  data  falls  below  a  predetermined 
threshold  criteria. 

In  any  inversion,  a  reasonable  guess  of  the  velocity  structure  in  the  starting  model  is  very  important .  To  test  the 
appropriateness  of  our  3-D  a  priori  model,  we  relocated  a  set  of  44  Calibration  Event  Bulletin  events  using  arrival 
times  form  14  regional  stations.  In  this  analysis,  our  3-D  model  resulted  in  smaller  rms  residuals  for  66  percent 
of  the  events.  Since  ground  truth  locations  are  not  available  for  these  events  and  the  depths  are  poorly  constrained 
by  the  regional  data,  further  analysis  was  performed  on  a  larger  data  set.  In  an  extension  of  this  study,  we 
relocated  200  events  from  the  Engdahl  el  al.  (1998)  database  using  our  current  3-D  model  (Fiugre  5). 


Figure  5.  A  set  of  200  events  from  the  Engdahl  et  al.  (1998)  database  used  in  an  preliminary  test  of  our  tomographic 
inversion  plan. 

We  chose  events  that  were  reported  in  the  catalog  to  have  standard  error  in  position  of  less  than  or  equal  to  5  km, 
as  determined  from  the  diagonal  elements  of  the  covariance  matrix  (Engdahl  el  al,  1998).  The  event  list  was 
further  narrowed  to  optimize  spatial  distribution  and  to  exclude  deep  events  (maximum  depth  was  163  km).  A 
subset  of  thirty  regional  stations,  also  selected  for  optimal  spatial  distribution  as  well  as  for  quantity  of  data 
reported  in  the  bulletin  for  each  station,  was  used  to  locate  each  event.  We  relocated  each  event  with  the  Weston 


3-D  velocity  model  of  Pakistan  and  then  with  the  1-D  IASP91  model.  This  test  was  performed  for  three  different 
cases  relative  to  the  Engdahl  et  al.  (1998)  hypocenter  solution:  (1)  constrained  to  a  10  km  radius  in  X-Y  space 
and  +/-  20  km  in  Z,  (2)  a  fixed  depth  constrained  to  +/-  5  km,  but  free  in  X-Y  space,  and  (3)  a  free  solution  in  all 
directions. 

The  hypocenter  solution  for  the  3-D  model  has  a  smaller  mis  residual  for  approximately  75%  of  the  events  in 
both  variations  of  the  constrained  cases  (1  and  2)  and  it  has  smaller  residuals  for  just  over  60%  of  the  events  in  the 
free  solution  (3)  (see  Table  3).  Average  standard  error  reductions  for  the  set  of  events  are  also  shown  to  be  1 30— 
170%  greater  than  average  standard  error  increases  for  the  3-D  hypocenter  solutions  compared  with  the  1-D 
solutions.  These  results  show  that  the  3-D  model  is  generally  fitting  the  data  better  than  the  1-D  model.  A 
smaller  rms  error  can  be  an  indication  of  a  better  location;  however,  3-D  structures  may  introduce  biases  that  lead 
to  smaller  rms  errors.  Therefore,  small  rms  residuals  are  not  an  absolute  indicator  of  a  better  location.  Other  tests 
(such  as  jackknifing,  as  discussed  in  the  previous  section)  are  required  to  demonstrate  that  locations  calculated 
using  the  3-D  model  are  more  accurate  than  the  1-D  solutions. 


Location 

Constraints 

No.  of  3-D  Events 
with  smaller  rms 
residuals  than  1-D 

No.  of  3-D  Events 
with  larger  rms 
residuals  than  1-D 

Average 
Standard  Error 
Reduction 

Average 
Standard  Error 
Increase 

Case  1:  XYZ* 

148 

51 

0.4042 

0.2098 

Case  2:  Z 

148 

52 

0.3920 

0.2298 

Case  3:  None 

121 

79 

0.2766 

0.2118 

*One  best  hypocenter  solution  was  found  outside  our  model  region  and  was  removed  from  this  evaluation 


Table  3.  Comparison  of  rms  residuals  for  preliminary  tomographic  analysis  of 200  events  from  the  Engdahl,  et  al. 
(1998)  database.  These  events  were  relocated  using  both  the  3-D  regional  model  and  the  1-D  reference  IASP91 
model.  Three  different  hypocenter  contraint  conditions  were  imposed  on  the  inversion. 

To  obtain  a  preview  of  the  magnitude  and  spatial  distribution  of  travel  time  residuals  over  the  model  space,  we 
plotted  the  residuals  for  each  station-event  pair  at  the  great-circle  midpoint  of  the  ray.  This  is  only  an  informal 
spatial  approximation,  as  the  residuals  are  distributed  over  the  entire  raypath.  Therefore,  interpretations  and 
conclusions  drawn  from  this  representation  of  the  results  are  limited  and  must  take  the  nature  of  the  midpoint 
approximation  into  consideration.  We  can  see  from  Figures  6  and  7  that  locations  determined  using  the  3-D 
model  result  in  overall  smaller  residuals  than  those  determined  by  the  IASP91  model.  In  both  cases,  the 
magnitudes  and  sign  of  the  residuals  appear  to  exhibit  a  spatial  correlation  across  some  regions.  For  most  regions, 
it  appears  that  the  3-D  model  already  does  a  better  job  of  approximating  the  true  velocity  structure  than  IASP91, 
as  expected.  However,  we  anticipate  that  full  implementation  of  the  tomographic  inversion,  which  will  update  the 
a  priori  model  along  the  full  raypath  between  each  station  and  event,  will  improve  the  regions  that  are  less  well 
constrained  in  the  preliminary  model.  This  initial  test  demonstrates  the  functionality  of  the  majority  of 
components  required  for  tomographic  inversion.  It  also  produces  a  "preview"  of  the  spatial  distribution  of  travel 
time  residuals  over  the  model  space. 


-3  00  -2  25  -1  50  -0  75  0  00  0.75  1  50  2  25  3  00 
Residual  <s) 


Figure  6.  Residuals  of  200  events  relocated  using  the  Weston  3-D  velocity  model  of  Pakistan.  The  solutions  were 
constrained  to  within  a  10  km  radius  in  XY  space  and  to  +/-  20  km  in  depth  of  the  Engdahl  el  al.  (1998)  solutions. 
Residuals  are  plotted  at  the  great  circle  midpoint  of  the  raypath  and  are  averaged  for  points  located  within  a  9  km 
grid  spacing.  Areas  not  resolved  from  this  small  set  of  the  data  are  white.  Large  residuals  (either  positive  or 
negative)  are  indicated  in  dark  shades,  whereas  light  shades  indicate  residuals  approaching  zero.  Highs  and  lows 
are  indicated  by  either  a"+"  or  sign.. 


•3.00  -2.25  -1.50  -0.75  0  00  0.75  1.50  2.26  3.03 


Residual  (s) 


Figure  7.  Residuals  from  the  same  set  of  events  shown  in  Figure  6  relocated  using  the  IASP91  1-D  travel  time  tables. 


CONCLUSIONS  AND  RECOMMENDATIONS 


During  the  past  year  we  have  made  significant  progress  towards  producing  a  fully  3-D  tomographic  velocity 
model  of  the  Pakistan  region  and  surrounding  areas.  We  have  completed  the  development  and  testing  of  a  3-D 
starting  velocity  model  and  have  established  a  database  of  seismic  events  that  we  are  using  as  the  input  travel 
times  to  the  inversion.  We  continue  to  evaluate  the  quality  of  each  of  our  input  data  sets  to  ensure  the  highest 
quality  data  is  available  for  our  inversion.  In  this  paper  we  have  described  the  functionality  and  demonstrated  the 
application  of  key  computational  modules  necessary  to  support  the  final  tomographic  inversion.  In  each  case,  we 
have  completed  various  performance  tests  on  the  computational  accuracy. 

Work  to  be  completed  during  the  final  year  of  this  contract  includes  the  final  development  and  implementation  of 
the  joint  inversion  technique  and  establishment  of  appropriate  convergence  criteria.  Once  these  modules  are 
completed  and  tested,  we  will  perform  a  full  tomographic  inversion  to  produce  an  accurate  regional  3-D  velocity 
model  for  the  Pakistan  region  and  the  consequent  improved  hypocenter  estimates. 
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